﻿using NumericalMatrix.Interfaces;
using NumericalMatrix.Exceptions;
using NumericalMatrix.GenericInterfaces;
using System;

namespace NumericalMatrix.Services
{
    //, MatrixType
    public class LUBlock<ElementType> : ILUBlock<ElementType> //where ElementType : IMatrix<ElementType>
    {
        #region Fields

        IMatrix<ElementType> l;
        IMatrix<ElementType> u;

        IMatrix<ElementType> matrix;

        IGroupOperations<ElementType> group;
        private Matrix matrix_2;

        #endregion

        #region Constructors

        public LUBlock(IMatrix<ElementType> matrix)
        {
            if (matrix.Height != matrix.Width)
            {
                //TODO: create MatrixWrongDimension exception
                throw new MatrixException("wrong dimension");
            }

            this.matrix = matrix;
        }

        public LUBlock(Matrix matrix_2)
        {
            // TODO: Complete member initialization
            this.matrix_2 = matrix_2;
        }

        #endregion

        #region interface members

        #region ILUBlock members

        public void CleanData()
        {
            // TODO make Matrix disposable
            L = null;
            U = null;
        }

        #endregion

        #region ILUInvertible members

        public IMatrix<ElementType> InvertLU()
        {
            throw new Exception();
            //var width = matrix.Width;
            //var height = matrix.Height;

            //IMatrix<ElementType> S;

            //if (typeof(ElementType) == typeof(double))
            //{
            //    S = matrix.CreateDoubleInstance(height, width) as IMatrix<ElementType>;
            //}
            //if (typeof(ElementType) == typeof(GroupOperators))
            //{
            //    S = matrix.CreateInstance(height, width) as IMatrix<ElementType>;
            //}

            
            //IMatrix<ElementType> x;
            //IMatrix<ElementType> y;

            //for (int i = 0; i < width; i++)
            //{
            //    Matrix2 b = new Matrix2(height, 1);
            //    b[i, 0] = 1;

            //    // TODO: generilize this extension methods

            //    y = L.LinSolveLtriang(b);
            //    x = U.LinSolveRtriang(y);
                
                
            //    //y.WriteToFile("y.txt");
            //    //(RH * x).WriteToFile("RH_X.txt");

            //    for (int j = 0; j < width; j++)
            //    {
            //        S[j, i] = x[j, 0];
            //    }

            //}

            //return S;
        }

        #endregion

        #region ILUFactorisable members

        public IMatrix<ElementType> L
        {
            get;
            private set;
        }

        public IMatrix<ElementType> U
        {
            get;
            private set;
        }

        public void LUFactorise()
        {
            var width = matrix.Width;
            var height = matrix.Height;

            L = matrix.CreateInstance(height, width) as IMatrix<ElementType>;
            U = matrix.CreateInstance(height, width) as IMatrix<ElementType>;

            for (int i = 0; i < width; i++)
                L[i, i] = matrix.GroupOperations.NeutralElement;
            
            //iteration by L rows
            for (int i = 0; i < width; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    ElementType s = group.ZeroElement;
                    for (int k = 0; k < j; k++)
                    {
                        var multipl = group.Multiplication(L[i, k], U[k, j]);
                        s = group.Addition(s, multipl);
                        //s += group.Multiplication(L[i, k], U[k, j]);
                    }

                    var subtract = group.Subtraction(matrix[i, j], s);
                    L[i, j] = group.Division(subtract, U[j, j]);
                    //L[i, j] = (matrix[i, j] - s) / U[j, j];
                }

                for (int j = i; j < width; j++)
                {
                    ElementType s = group.NeutralElement;
                    for (int k = 0; k < i; k++)
                    {
                        var multipl = group.Multiplication(L[i, k], U[k, j]);
                        s = group.Addition(s, multipl);
                        //s += L[i, k] * U[k, j];
                    }

                    U[i, j] = group.Subtraction(matrix[i, j], s);
                    //U[i, j] = matrix[i, j] - s;
                }
            }
        }

        #endregion

        #endregion
    }
}
